from concurrent.futures import ThreadPoolExecutor
from concurrent.futures import ProcessPoolExecutor
try:
    import gdal
except:
    from osgeo import gdal
try:
    import gdalconst
except:
    from osgeo import gdalconst
try:
    import osr
except:
    from osgeo import osr
from xml.dom import minidom
import pandas as pd
import numpy as np
import argparse
import datetime
import queue
from tqdm import tqdm
import glob
import os
import re

os.environ['PROJ_LIB'] = r"/home/fengyangyang/python_env_1/share/proj"

class BoundedThreadPoolExecutor(ThreadPoolExecutor):
    def __init__(self, max_workers=None, thread_name_prefix=''):
        super().__init__(max_workers,thread_name_prefix)
        self._work_queue = queue.Queue(max_workers * 2)

# 进度条
def progress_cb(complete, message, cb_data):
    '''Emit progress report in numbers for 10% intervals and dots for 3%'''
    if int(complete*100) % 10 == 0:
        print(f'{complete*100:.0f}', end='', flush=True)
    elif int(complete*100) % 3 == 0:
        print(f'{cb_data}', end='', flush=True)

# 获得编码
class get_gy100_code_by_rule:

	# （首位编号为0）获得图幅编号第1号位字母
	first_letters = ["Q", "R", "S", "T", "U"]
	# 中国区度带号列表
	zones = [f"{zone}" for zone in range(43, 54)]
	point = (500000, 1900000)

	def __init__(self, mgrs):
		self.mgrs = mgrs
		self.zone = mgrs[:2]
		self.letter_1bit = mgrs[2]
		self.letter_2bit = mgrs[3]
		self.letter_3bit = mgrs[4]

	# 获得图幅编号第3号位字母列表
	@staticmethod
	def get_y_letters():
		yletters = [chr(i) for i in range(65, 91)]
		yletters.remove("I")
		yletters.remove("O")
		yletters.remove("W")
		yletters.remove("X")
		yletters.remove("Y")
		yletters.remove("Z")
		yletter_count = len(yletters)
		letter_ys = [yletters[i % yletter_count] for i in range(100)]
		three1_letters = letter_ys[letter_ys.index("T"):]
		three0_letters = letter_ys[letter_ys.index("C"):]
		return [three0_letters, three1_letters]

	def get_zone_letter2bits(self):
		""""""
		second_letters = self.get_x_letters()
		x_letters = []
		for zone in self.zones:
			for index in range(8):
				letter_2bit = second_letters.pop(0)
				x_letters.append([zone, letter_2bit])
		return x_letters

	@classmethod
	def get_letter13bits(cls):
		""""""
		[letter3bits_0, letter3bits_1] = cls.get_y_letters()
		letters_13bits_0 = [(None, None)]
		for first_letter in cls.first_letters:
			for index in range(10):
				if index == 0 and first_letter != "Q":
					letters_13bits_0[-1] = (letters_13bits_0[-1], [first_letter, letters_13bits_0[-1][1]])
				else:
					letters_13bits_0.append([first_letter, letter3bits_0.pop(0)])

		letters_13bits_1 = [(None, None)]
		for first_letter in cls.first_letters:
			for index in range(10):
				if index == 0 and first_letter != "Q":
					letters_13bits_1[-1] = (letters_13bits_1[-1], [first_letter, letters_13bits_1[-1][1]])
				else:
					letters_13bits_1.append([first_letter, letter3bits_1.pop(0)])

		letters_13bitss = [letters_13bits_0, letters_13bits_1]
		return letters_13bitss

	def get_x_index(self):
		""""""
		x_letters = self.get_zone_letter2bits()
		x_index = x_letters.index([self.zone, self.letter_2bit])
		return x_index

	def get_y_index(self):
		letter13bits = self.get_letter13bits()[int(self.zone) % 2]
		y_index = None
		for index, letters_13bit in enumerate(letter13bits):
			if isinstance(letters_13bit, list):
				if [self.letter_1bit, self.letter_3bit] == letters_13bit:
					y_index = index
					break
			else:
				if [self.letter_1bit, self.letter_3bit] in letters_13bit:
					y_index = index
					break
		return y_index

	# 获得图幅编号第2号位字母列表
	@staticmethod
	def get_x_letters():
		xletters = [chr(i) for i in range(65, 91)]
		xletters.remove("I")
		xletters.remove("O")
		xletter_count = len(xletters)
		letter_xs = [xletters[i % xletter_count] for i in range(100)]
		second_letters = letter_xs[letter_xs.index("A"):]
		return second_letters

	@property
	def gy100_code(self):
		""""""
		(x, y) = self.get_gy100_point()
		x_point_str = "0"*(2-len(f"{int(x//100000)}")) + f"{int(x//100000)}"
		y_point_str = "0"*(2-len(f"{int(y//100000)}")) + f"{int(y//100000)}"
		y_index = self.get_y_index()
		x_index = self.get_x_index()
		y_indexstr = "0" * (3-len(f"{y_index}")) + f"{y_index}"
		x_indexstr = "0" * (3-len(f"{x_index}")) + f"{x_index}"
		gy100_code = f"{self.mgrs}_{y_indexstr}_{x_indexstr}_100_{x_point_str}_{y_point_str}"
		return gy100_code

	# 根据mgr获得100km的point
	def get_gy100_point(self):
		""""""
		# 格投影度带参考
		base_y_index = 2
		base_x_indexs = [4 + 8 * i for i in range(len(self.zones))]
		zone_letter2bits = self.get_zone_letter2bits()
		for x_index in base_x_indexs:
			[zone, letter2bit] = zone_letter2bits[x_index]
			if zone == self.zone:
				base_x_index = x_index
		y_index = self.get_y_index()
		x_index = self.get_x_index()
		x = self.point[0] + (x_index - base_x_index) * 100000
		y = self.point[1] + (y_index - base_y_index) * 100000
		return (x, y)

# 依据Sentinel2格网Code计算GY100格网的Code
def get_gy100_code(tile):
	return get_gy100_code_by_rule(tile).gy100_code

# gy100km转gy10km
def get_gy10_codes(gy100):
	""""""
	gy100_params = gy100.split("_")
	gy10_codes = []
	for x in range(10):
		for y in range(10):
			gy10_params = [
				f"{gy100_params[0]}",
				f"{gy100_params[1]}",
				f"{gy100_params[2]}",
				f"{y}",
				f"{x}",
				f"10",
				f"{gy100_params[-2]}{x}",
				f"{int(gy100_params[-1]) * 10 - y}",
			]
			gy10km = "_".join(gy10_params)
			gy10_codes.append(gy10km)
	return gy10_codes

# 根据gy100_code获得左上角坐标点
def get_gy100_point(tile):
	gy100_code = get_gy100_code(tile)
	gy100_params = gy100_code.split("_")
	x = int(gy100_params[-2]) * 100000
	y = int(gy100_params[-1]) * 100000
	gy100_point = [x, y]
	return gy100_point

# 读取xml文件中的数据
def read_xml_data(file_path, select_node):
	""""""
	domTree = minidom.parse(file_path)
	collection = domTree.documentElement
	select_element = [collection]
	for node_dict in select_node:
		tag_name = node_dict["tag_name"]
		if node_dict["attribute"] is None:
			select_element = [
				element_sub
				for element in select_element
			    for element_sub in element.getElementsByTagName(tag_name)
			]
		else:
			attribute_name = node_dict["attribute"]["name"]
			attribute_value = node_dict["attribute"]["value"]
			select_element = [
				element_sub
				for element in select_element
				for element_sub in element.getElementsByTagName(tag_name)
				if element_sub.getAttribute(attribute_name) == attribute_value
			]
	datas = [
		element.firstChild.data
		for element in select_element
	]
	return datas

# 生成指定编号的gy100栅格文件路径
def get_gy100km_raster(
    gykm_dir, tile, year, month, period, band, sentinel_dir=None):
	"""
	:param gykm_dir: gykm 数据路径
	:param tile: "51TVN"
	:param year: 2021
	:param month: 7
	:param period: 1501
	:param band: 5
	:param sentinel_dir: /83To82/all_data/he_nan/2021_12/3A/202112/3000/50RLV/SENTINEL2X_20211230-000000-000_L3A_T50RLV_C_V1-2
	:return:
	"""
	# 获得gy100目录
	gy100 = get_gy100_code(tile)
	date = datetime.datetime(year, month, 1)
	yearmonth = date.strftime("%Y%m")
	gy100km_dir = os.path.join(
		gykm_dir, f"326{tile[:2]}", f"{tile[:3]}", f"100KM", f"{gy100}",
		f"Image_Data", f"{yearmonth}", f"L3A_TIFF", f"{yearmonth}_{period}",
	)
	try:
		if not os.path.exists(gy100km_dir):
			os.makedirs(gy100km_dir)
	except:
		pass
	# 获得gy100 raster 路径
	if sentinel_dir is not None:
		cloud = read_cloud(sentinel_dir)
		gy100km_raster_name = f"{gy100}_{yearmonth}_{period}_{cloud}_B{band}.tif"
		gy100km_raster_path = os.path.join(gy100km_dir, gy100km_raster_name)
	else:
		gy100km_raster_namere = f"{gy100}_{yearmonth}_{period}_*_B{band}.tif"
		gy100km_raster_pathre = os.path.join(gy100km_dir, gy100km_raster_namere)
		gy100km_raster_paths = glob.glob(gy100km_raster_pathre)
		if len(gy100km_raster_paths) == 1:
			gy100km_raster_path = gy100km_raster_paths[0]
		else:
			gy100km_raster_path = None
	return gy100km_raster_path

# 获得含云量
def read_cloud(sentinel_dir):
	""""""
	# print(sentinel_dir)
	sentinel_name = os.path.basename(sentinel_dir)
	l1c_namere = f"S2[AB]_MSIL1C_.*\.SAFE"
	l2a_namere = f"SENTINEL2[AB]_.*_L2A_T.*_C_V1-0"
	l3a_namere = f"SENTINEL2[ABX]_.*_L3A_T.*_C_V1-2"
	if re.match(l1c_namere, sentinel_name):
		xml_name = f"MTD_MSIL1C.xml"
		xml_path = os.path.join(sentinel_dir, xml_name)
		select_node = [
			{
				"tag_name": "n1:Quality_Indicators_Info",
				"attribute": None
			},
			{
				"tag_name": "Cloud_Coverage_Assessment",
				"attribute": None
			}
		]
		# 读取xml文件中的数据
		datas = read_xml_data(xml_path, select_node)
		cloud = f"{round(float(datas[0]))}"
		cloud = f"0{cloud}" if len(cloud) == 1 else cloud
	elif re.match(l2a_namere, sentinel_name):
		xml_name = f"{sentinel_name}_MTD_ALL.xml"
		xml_path = os.path.join(sentinel_dir, xml_name)
		select_node = [
			{
				"tag_name": "Current_Product",
				"attribute": None
			},
			{
				"tag_name": "Product_Quality_List",
				"attribute": {"name": "level", "value": "N3"}
			},
			{
				"tag_name": "Global_Index_List",
				"attribute": None
			},
			{
				"tag_name": "QUALITY_INDEX",
				"attribute": {"name": "name", "value": "CloudPercent"}
			}
		]
		# 读取xml文件中的数据
		# datas = read_xml_data(xml_path, select_node)
		# cloud = f"{round(float(datas[0]))}"
		# cloud = f"0{cloud}" if len(cloud) == 1 else cloud
		cloud = "00"
	elif re.match(l3a_namere, sentinel_name):
		xml_name = f"{sentinel_name}_MTD_ALL.xml"
		xml_path = os.path.join(sentinel_dir, xml_name)
		select_node = [
			{
				"tag_name": "Current_Product",
				"attribute": None
			},
			{
				"tag_name": "Product_Quality_List",
				"attribute": {"name": "level", "value": "N3"}
			},
			{
				"tag_name": "Global_Index_List",
				"attribute": None
			},
			{
				"tag_name": "QUALITY_INDEX",
				"attribute": {"name": "name", "value": "CloudPercent"}
			}
		]
		# 读取xml文件中的数据
		# datas = read_xml_data(xml_path, select_node)
		cloud = "00"

	# 读取xml文件中的数据
	# datas = read_xml_data(xml_path, select_node)
	# cloud = f"{round(float(datas[0]))}"
	# cloud = f"0{cloud}" if len(cloud) == 1 else cloud
	return cloud

	""""""
	sentinel_name = os.path.basename(sentinel_dir)
	l1c_namere = f"S2[AB]_MSIL1C_.*\.SAFE"
	l2a_namere = f"SENTINEL2[AB]_.*_L2A_T.*_C_V1-0"
	l3a_namere = f"SENTINEL2[ABX]_.*_L3A_T.*_C_V1-2"
	if f"{band}" == f"99":
		if re.match(l1c_namere, sentinel_name):
			flag = "1C_MASK"
		elif re.match(l2a_namere, sentinel_name):
			flag = "2A_MASK"
		elif re.match(l3a_namere, sentinel_name):
			flag = "3A_MASK"
		else:
			flag = False
	else:
		flag = True
	return flag

# 获得gy10km切片参数
def get_gy10_points(tile):
    gy100_code = get_gy100_code(tile)
    gy10_codes = get_gy10_codes(gy100_code)
    gy10_points = []
    for gy10_code in gy10_codes:
        gy10_infos = gy10_code.split("_")
        left = int(gy10_infos[-2]) * 10000.0
        top = int(gy10_infos[-1]) * 10000.0
        point = [left, top]
        gy10_points.append(point)
    return gy10_points

# delete old gy100km data 
def delete_old_gy100km(
	tiles, bands, year, month, period, gykm_dir):
	yearmonth = f"{year}" + "0"*(len(f"{month}")<2) + f"{month}"
	for index in tqdm(range(len(tiles))):
		tile = tiles[index]
		gy100_code = get_gy100_code(tile)
		gy100_dir = os.path.join(
			gykm_dir, 
			f"326{tile[:2]}", 
			f"{tile[:3]}", 
			f"100KM", 
			f"{gy100_code}",
			f"Image_Data", 
			f"{yearmonth}", 
			f"L3A_TIFF", 
			f"{yearmonth}_{period}")
		if not os.path.exists(gy100_dir):
			continue
		for band in bands:
			gy100_namere_params = [
				gy100_code, 
				yearmonth , 
				f"{period}", 
				"*", 
				f"B{band}.tif"
			]
			gy100_namere = "_".join(gy100_namere_params)
			gy100_pathre = os.path.join(gy100_dir, gy100_namere)
			for path in glob.glob(gy100_pathre):
				os.remove(path)
				# print(f"{datetime.datetime.now()}>删除:{path}")
	return

# delete old gy10km data 
def delete_old_gy10km(
	tiles, bands, year, month, period, gykm_dir):
	yearmonth = f"{year}" + "0"*(len(f"{month}")<2) + f"{month}"
	for index in tqdm(range(len(tiles))):
		tile = tiles[index]
		gy100_code = get_gy100_code(tile)
		gy10_codes = get_gy10_codes(gy100_code)
		for gy10_code in gy10_codes:
			gy10_dir = os.path.join(
				gykm_dir, 
				f"326{tile[:2]}", 
				f"{tile[:3]}", 
				f"10KM", 
				f"{gy100_code}",
				f"{gy10_code}",
				f"Image_Data", 
				f"{yearmonth}", 
				f"L3A_TIFF", 
				f"{yearmonth}_{period}")
			if not os.path.exists(gy10_dir):
				continue
			for band in bands:
				gy10_namere_params = [
					gy10_code, 
					yearmonth , 
					f"{period}", 
					"*", 
					f"B{band}.tif"
				]
				gy10_namere = "_".join(gy10_namere_params)
				gy10_pathre = os.path.join(gy10_dir, gy10_namere)
				for path in glob.glob(gy10_pathre):
					os.remove(path)
					# print(f"{datetime.datetime.now()}>删除:{path}")
	return

# 统计合成3a所需那些日期的2a
def get_l2a_info(tiles, folder):
	l1c_namere = f"S2[AB]_MSIL1C_.*\.SAFE"
	l2a_namere = f"SENTINEL2[AB]_.*_L2A_T.*_C_V1-0"
	l3a_namere = f"SENTINEL2[ABX]_.*_L3A_T.*_C_V1-2"
	results = []
	for index in tqdm(range(len(tiles))):
		tile = tiles[index]
		tile_folder = os.path.join(folder, tile)
		if not os.path.exists(tile_folder):
			# print(f"{tile}:无数据，请合成或补充数据！")
			results.append(f"{tile}:无数据，请合成或补充数据！")
			continue
		flag = True
		for name in os.listdir(tile_folder):
			if re.match(l1c_namere, name):
				sentinel_dir = os.path.join(tile_folder, name)
				date = name.split("_")[2][:8]
				# print(f"{tile}:type(L1C),{date}")
				results.append(f"{tile}:type(L1C),{date}")
				flag = False
				break
			elif re.match(l2a_namere, name):
				sentinel_dir = os.path.join(tile_folder, name)
				date = name.split("_")[1][:8]
				# print(f"{tile}:type(L2A),{date}")
				results.append(f"{tile}:type(L2A),{date}")
				flag = False
				break
			elif re.match(l3a_namere, name):
				sentinel_dir = os.path.join(tile_folder, name)
				xml_name = f"{name}_MTD_ALL.xml"
				xml_path = os.path.join(sentinel_dir, xml_name)
				select_node = [
					{
						"tag_name": "Contributing_Products_List",
						"attribute": None,
					},
					{
						"tag_name": "Product",
						"attribute": None,
					},
					{
						"tag_name": "PRODUCT_ID",
						"attribute": None,
					},
				]
				l2a_names = read_xml_data(xml_path, select_node)
				dates = [
					l2a_name.split("_")[1][:8]
					for l2a_name in l2a_names]
				# print(f"{tile}:type(L3A),{dates}")
				results.append(f"{tile}:type(L3A),{dates}")
				flag = False
				break
		if flag:
			# print(f"{tile}:无数据，请合成或补充数据！")
			results.append(f"{tile}:无数据，请合成或补充数据！")
			pass
	return results

# create l1c mask raster
def get_l1c_mask_raster(raster_path):
    name = os.path.basename(raster_path)
    tile = name.split("_")[5][1:]
    spatial = osr.SpatialReference()
    spatial.ImportFromEPSG(int(f"326{tile[:2]}"))
    xsize, ysize = 10000, 10000
    array = np.zeros((xsize, ysize), dtype=np.uint16)
    [left, top] = get_gy100_point(tile)
    geotransform = (left, 10.0, 0.0, top, 0.0, -10.0)
    gdal.AllRegister()
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(
        raster_path,
        xsize=xsize,
        ysize=ysize,
        bands=1, 
        eType=gdalconst.GDT_UInt16
    )
    dataset.SetSpatialRef(spatial)
    dataset.SetGeoTransform(geotransform)
    band = dataset.GetRasterBand(1)
    band.WriteArray(array)
    return

# get 1c raster path
def get_l1c_raster_path(sentinel_dir, band):
    folder_name = os.path.basename(sentinel_dir)
    tile = folder_name.split("_")[5][1:]
    date = folder_name.split("_")[2]
    if f"{band}" == f"99":
        raster_name = folder_name.replace(".SAFE", "_MASK.tif")
        raster_path = os.path.join(sentinel_dir, raster_name)
        if not os.path.exists(raster_path):
            get_l1c_mask_raster(raster_path)
    else:
        band = "0"*(len(f"{band}")==1) + f"{band}"
        raster_name = f"T{tile}_{date}_B{band}.jp2"
        raster_pathre = os.path.join(
            sentinel_dir, 
            "GRANULE",
            f"L1C_T{tile}_A*_{date[:8]}T*",
            "IMG_DATA", 
            raster_name)
        raster_path = glob.glob(raster_pathre)[0]
    return raster_path 

# get 2a raster path
def get_l2a_raster_path(sentinel_dir, band):
    sentinel_name = os.path.basename(sentinel_dir)
    if f"{band}" == f"99":
        mask_name = f"{sentinel_name}_MG2_R1.tif"
        sentinel_path = os.path.join(
            sentinel_dir, "MASKS", mask_name)
    else:
        band_name = f"{sentinel_name}_FRE_B{band}.tif"
        sentinel_path = os.path.join(sentinel_dir, band_name)
    return sentinel_path

# get 3a raster path
def get_l3a_raster_path(sentinel_dir, band):
    sentinel_name = os.path.basename(sentinel_dir)
    if f"{band}" == f"99":
        mask_name = f"{sentinel_name}_FLG_R1.tif"
        sentinel_path = os.path.join(
            sentinel_dir, "MASKS", mask_name)
    else:
        band_name = f"{sentinel_name}_FRC_B{band}.tif"
        sentinel_path = os.path.join(sentinel_dir, band_name)
    return sentinel_path

# get gy10 raster paths
def get_gy10_raster_paths(
    tile, band, year, month, period, gykm_dir):
    gy100_code = get_gy100_code(tile)
    yearmonth = f"{year}" + "0"*(len(f"{month}")<2) + f"{month}"
    gy10_codes = get_gy10_codes(gy100_code)
    gy10_paths = []
    for gy10_code in gy10_codes:
        gy10_dir = os.path.join(
            gykm_dir, 
            f"326{tile[:2]}", 
            f"{tile[:3]}", 
            f"10KM", 
            f"{gy100_code}",
            f"{gy10_code}",
            f"Image_Data", 
            f"{yearmonth}", 
            f"L3A_TIFF", 
            f"{yearmonth}_{period}")
        if not os.path.exists(gy10_dir):
            os.makedirs(gy10_dir)
        gy10_name = f"{gy10_code}_{yearmonth}_{period}_00_B{band}.tif"
        gy10_path = os.path.join(gy10_dir, gy10_name)
        gy10_paths.append(gy10_path)
    return gy10_paths

# get gy100 raster paths
def get_gy100_raster_path(
    tile, band, year, month, period, gykm_dir):
    gy100_code = get_gy100_code(tile)
    yearmonth = f"{year}" + "0"*(len(f"{month}")<2) + f"{month}"
    gy100_dir = os.path.join(
        gykm_dir, 
        f"326{tile[:2]}", 
        f"{tile[:3]}", 
        f"100KM", 
        f"{gy100_code}",
        f"Image_Data", 
        f"{yearmonth}", 
        f"L3A_TIFF", 
        f"{yearmonth}_{period}")
    if not os.path.exists(gy100_dir):
        os.makedirs(gy100_dir)
    gy100_name = f"{gy100_code}_{yearmonth}_{period}_00_B{band}.tif"
    gy100_path = os.path.join(gy100_dir, gy100_name)
    return gy100_path

# get band raster path
def get_band_raster_path(tile, band, folder):
	l1cnamere = "^S2[AB]_MSIL1C_.*\.SAFE$"
	l2anamere = "^SENTINEL2[AB]_.*_L2A_T.*_C_V1-0$"
	l3anamere = "^SENTINEL2[ABX]_.*_L3A_T.*_C_V1-2$"
	raster_path = None
	tile_folder = os.path.join(folder, tile)
	if not os.path.exists(tile_folder):
		return raster_path
	# flag = True
	for name in os.listdir(tile_folder):
		if re.match(l1cnamere, name):
			sentinel_dir = os.path.join(tile_folder, name)
			raster_path = get_l1c_raster_path(sentinel_dir, band)
			# flag = False
			break
		elif re.match(l2anamere, name):
			sentinel_dir = os.path.join(tile_folder, name)
			raster_path = get_l2a_raster_path(sentinel_dir, band)
			# flag = False
			break
		elif re.match(l3anamere, name):
			sentinel_dir = os.path.join(tile_folder, name)
			raster_path = get_l3a_raster_path(sentinel_dir, band)
			# flag = False
			break
	# if flag:
	# 	print(f"{datetime.datetime.now()}>缺少数据:{tile}, {band}")
	return raster_path

# get all gykm raster files
def get_gykm_raster(
	tiles, bands, year, month, period, folder, gykm_dir):
	gy100_args = []
	gy10_args = []
	q = queue.Queue(maxsize=1)
	q.put(0)
	for tile in tiles:
		for band in bands:
			band_path = get_band_raster_path(tile, band, folder)
			point = get_gy100_point(tile)
			gy100_path = get_gy100_raster_path(
				tile, band, year, month, period, gykm_dir)
			gy100_args.append([band_path, point, gy100_path, q])
			points = get_gy10_points(tile)
			gy10_paths = get_gy10_raster_paths(
				tile, band, year, month, period, gykm_dir)
			gy10_args.append([gy100_path, points, gy10_paths, q])       
	# with ThreadPoolExecutor(max_workers=50) as executor:
	with BoundedThreadPoolExecutor(max_workers=50) as executor:
		results = executor.map(get_gy100_raster, gy100_args)   
	gy100_errors = [ 
		tile
		for result, tile in zip(results, tiles) 
		if not result] 
	# with ThreadPoolExecutor(max_workers=50) as executor:
	with BoundedThreadPoolExecutor(max_workers=50) as executor:
		results = executor.map(get_gy10_rasters, gy10_args)
	gy10_errors = [ 
		tile
		for result, tile in zip(results, tiles) 
		if not result] 
	return gy100_errors, gy10_errors

# get gy100km raster file
def get_gy100_raster(args):
	[source_path, point, target_path] = args
	if source_path is None:
		name = os.path.basename(target_path)
		tile = name.split("_")[0]
		return tile
	[left, top] = point
	right = left + 10 * 10000
	bottom = top - 10 * 10000
	outputBounds = (left, bottom, right, top)
	gdal.Warp(
		target_path,
		source_path,
		format='GTiff', 
		xRes = 10,
		yRes = 10,
		outputBounds=outputBounds, 
		resampleAlg=gdalconst.GRA_NearestNeighbour,
		multithread=20,
		# callback=progress_cb,
		# callback_data='.'
	)
	return None

# get gy10km raster files
def get_gy10_rasters(args):
	[source_path, points, target_paths] = args
	if not os.path.exists(source_path):
		name = os.path.basename(source_path)
		tile = name.split("_")[0]
		return tile
	for point, target_path in zip(points, target_paths):
		[left, top] = point
		right = left + 10 * 1000
		bottom = top - 10 * 1000
		outputBounds = (left, bottom, right, top)
		gdal.Warp(
			target_path,
			source_path,
			format='GTiff', 
			xRes = 10,
			yRes = 10,
			outputBounds=outputBounds, 
			resampleAlg=gdalconst.GRA_NearestNeighbour,
			multithread=20,
			# callback=progress_cb,
			# callback_data='.'
		)
	return None

# start 
def main():
 
	parser = argparse.ArgumentParser(description='Process some args.')
	parser.add_argument(
		'--tiles', type=str, required=False, help="tiles", default=None)
	parser.add_argument(
		'--year', type=str, required=False, help="year", default=None)
	parser.add_argument(
		'--month', type=str, required=False, help="month", default=None)
	parser.add_argument(
		'--period', type=str, required=False, help="period", default=None)
	parser.add_argument(
		'--bands', type=str, required=False, help="bands", default=None)
	parser.add_argument(
		'--folder', type=str, required=False, help="folder", default=None)

	args = parser.parse_args()
	tiles = args.tiles
	year = args.year
	month = args.month
	period = args.period
	bands = args.bands
	folder = args.folder

	tiles = [tile.upper() for tile in tiles.split("-")] # 图幅号列表
	bands = [band for band in bands.split("-")]  # 波段号列表
	print(f"处理图幅列表({len(tiles)})：{tiles}")
	print(f"处理波段列表({len(bands)})：{bands}")

	"""只需配置一次即可"""
	gykm_dir = r"/data/jiabing/Automatic_interpretation"  # 正式
	# gykm_dir = r"/data/fengyangyang/down_test"  # 测试
	print("启动L2A数据检测结果！")
	results = get_l2a_info(tiles, folder)
	tiles = [
		result[:5]
		for result in results
		if not re.match(".*:无数据，请合成或补充数据！$", result)
	]
	print(f"有效图幅列表({len(tiles)}):{tiles}")
	print(f"删除GY100KM历史数据！")
	delete_old_gy100km(tiles, bands, year, month, period, gykm_dir)
	print(f"删除GY10KM历史数据！")
	delete_old_gy10km(tiles, bands, year, month, period, gykm_dir)
	print(f"启动GYKM切片任务！")
	
	print(f"准备切片参数：")
	gy100_args, gy10_args = [], []
	for tile in tiles:
		for band in bands:
			band_path = get_band_raster_path(tile, band, folder)
			point = get_gy100_point(tile)
			gy100_path = get_gy100_raster_path(
				tile, band, year, month, period, gykm_dir)
			gy100_args.append([band_path, point, gy100_path])
			points = get_gy10_points(tile)
			gy10_paths = get_gy10_raster_paths(
				tile, band, year, month, period, gykm_dir)
			gy10_args.append([gy100_path, points, gy10_paths]) 

	start_time = datetime.datetime.now()

	print(f"启动GY100KM切片：")
	with ProcessPoolExecutor(max_workers=50) as executor:
		executor.map(get_gy100_raster, gy100_args)  

	print(f"启动GY10KM切片：")
	with ProcessPoolExecutor(max_workers=50) as executor:
		executor.map(get_gy10_rasters, gy10_args)

	end_time = datetime.datetime.now()

	l2a_info = '\n'.join(results)
	print(f"2A数据检测结果：\n{l2a_info}")
	error_tiles = [
		result[:5]
		for result in results
		if re.match(".*:无数据，请合成或补充数据！$", result)
	]
	print(f"GYKM切片失败图幅:{error_tiles}")
	time_long = (end_time-start_time).seconds/60
	print(f"耗时：{time_long}分钟")
	print("END!")
	return


if __name__ == "__main__":
	# 3A数据切片程序入口
	main()
